Diffusive Counter Dispersion of Mass in Bubbly Media 
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We consider a liquid bearing gas bubbles in a porous medium. When gas bubbles are immovably 
trapped in a porous matrix by surface-tension forces, the dominant mechanism of transfer of gas mass 
becomes the diffusion of gas molecules through the liquid. Essentially, the gas solution is in local 
thermodynamic equilibrium with vapor phase all over the system, i.e., the solute concentration 
equals the solubility. When temperature and/or pressure gradients are applied, diffusion fluxes 
appear and these fluxes are faithfully determined by the temperature and pressure fields, not by the 
local solute concentration, which is enslaved by the former. We derive the equations governing such 
systems, accounting for thermodiffusion and gravitational segregation effects which are shown not 
to be neglected for geological systems — marine sediments, terrestrial aquifers, etc. The results are 
applied for the treatment of non-high-pressure systems and real geological systems bearing methane 
or carbon dioxide, where we find a potential possibility of the formation of gaseous horizons deep 
below a porous medium surface. The reported effects are of particular importance for natural 
methane hydrate deposits and the problem of burial of industrial production of carbon dioxide in 
deep aquifers. 

PACS numbers: 47.55.db, 66.10.C-, 92.40.Kf 



I. INTRODUCTION 

The diffusion of solute gases in liquids is a well-studied 
and relatively well- understood problem; see, e.g. 0, 
In the classical formulation of the problem, the diffusion 
flux of guest atoms or molecules is expressed in terms 
of various thermodynamic "forces" — concentration, pres- 
sure, and temperature gradients. Usually, the liquid-gas 
interface does not influence the nature of the bulk dif- 
fusion and determines only the boundary conditions for 
the flux. There exists, however, a vast class of systems, 
called "bubbly liquids", e.g. [H, where the liquid-gas in- 
terface maintains the gas-solution saturation all over the 
liquid volume and serves as a source/sink for the diffu- 
sion flux of the solute molecules (i.e., the diffusion flux 
does not change the solute concentration with time, but 
redistributes the mass between gas bubbles). Thus the 
liquid-gas interface plays an important role in the bulk 
diffusion of solutes. On the macroscopic scale, the diffu- 
sive mass transport of solutes may be very unusual and 
surprising in these systems. 

In the present study, we address bubbly liquids with 
immobile bubbles. These may be trapped by the sur- 
face forces in a porous matrix, therefore, for concrete- 
ness, we will consider liquids in porous media. Among 
numerous examples of such systems, which are of great 
practical importance, are the oil-bearing porous massifs, 
seabeds in organic carbon-rich marine sediments [IHH , 
aquifers 0, [Icj, peat-bogs and swamps. For these sys- 
tems, the presence of immobilized gas bubbles or liquid 
drops is well established experimentally. The respective 
solute gases include methane, carbon dioxide, oxygen, 
nitrogen, etc. Moreover, the problem of diffusion of the 
carbon dioxide in porous medium, filled with bubbly liq- 
uid (water) is directly related to the problem of burial 
of industrial release of CO2, e.g. [IJ 13 ■ I n Fig- [U we 



sketch a typical system, where the diffusion of a solute 
gas takes place in a bubbly liquid in a porous matrix. 

We will study diffusion transport on large space and 
time scales. Namely, we will be interested in the space 
scale much larger than the characteristic size of a gas 
bubble and in the time scale much larger than the char- 
acteristic relaxation time for gas dissolution in liquid; 
hence we assume that all over the volume, the solute 
gas is saturated in the solution. Indeed, for marine sed- 
iments the solution relaxation times measured by hours 
are small compared to the global evolution time scales, 
which can be as large as millions of years. To estimate 
the pore size I needed to trap a bubble of the compara- 
ble size we notice that the surface-tension forces, equal 
to (rl (a is the surface tension), should overwhelm the 
buoyancy force pu q gl 3 (pu q is the density of a liquid; g 
is the gravity). This yields the estimate, I < y/u/puqg, 
implying the maximal pore size I < 2.7 mm for water, 
which suggests trapping even for sands. We consider the 
immobilization of pore-size bubbles, because a big mov- 
ing bubble is unstable to splitting [lj]; it either stops or 
experiences splitting into smaller bubbles until their size 
becomes comparable to the pore size. We assume that 
our space scale is much larger than I. It is not straight- 
forward to estimate the relaxation time for the gas dis- 
solution; therefore, we just assume that the condition of 
the gas saturation in liquid is always fulfilled. 

Another important feature of bubbly liquids in porous 
media (relevant for many geological systems) is the pres- 
ence of a temperature gradient (usually due to the 
geothermal gradient) and a pressure gradient (due to 
the hydrostatic pressure). On the microscopic level, the 
temperature gradient causes the thermodiffusion (the so- 
called "Soret effect" |l|,|2|,|l4|,ll5|]), while the pressure gra- 
dient manifests the presence of the gravitational forces. 
Due to different mobilities and masses of the solute and 
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FIG. 1: Typical system of bubbly liquid filling the porous 
medium. The pressure, temperature, and parameters of the 
porous matrix may vary vertically, but are constant horizon- 
tally. The figure, with water above the bubble-rich sediments, 
sketches numerous geological systems. 



solvent molecules these may cause an additional mass 
flux. Although the importance of thermodiffusion for 
soil gas exchange has been reported in Ref . [l6[ , only the 
diffusion flux due to the concentration gradient (Fickian 
diffusion) has been investigated so far [^-Q. 

In the present study, we consider the diffusive trans- 
port of a solute gas in liquid in a porous medium with- 
out any significant through-flow in the system; that is, 
we consider only molecular and not hydrodynamic diffu- 
sion 0, [l(| Qjl [3 ■ We show that for large space and time 
scales, which assume the saturation condition for the so- 
lute gas, and under temperature and pressure gradients a 
novel and interesting phenomenon may be observed — the 
diffusive accumulation of gas in certain regions, caused 
by the non-Fickian transport of solutes. 

The paper is organized as follows. In Sec. [TT| the dif- 
fusion flux of a solute gas in bubbly liquid is considered 
in the presence of pressure and temperature gradients; 
the evolution equation for the gas concentration is ob- 
tained. In Sec. IIIII we demonstrate the existence of a 
novel phenomenon — the diffusive accumulation of solute 
gas in certain regions. We analyze this effect for real ge- 
ological systems — the methane or carbon dioxide bear- 
ing massifs in the presence of the geothermal gradient. 
Finally, in Sec. IIV1 we summarize our findings. Some 
technical details are given in the Appendixes. 



II. DIFFUSION TRANSPORT IN BUBBLY 
LIQUID WITH GRADIENTS OF 
THERMODYNAMIC PROPERTIES 

The diffusion current of guest particles in a host ma- 
terial in the presence of concentration and temperature 
gradients and an external potential generally has three 
components, each being the response to their respective 
driving force: The Fickian flux, due to the concentration 
gradient, the thermodiffusion flux, due to the tempera- 
ture gradient, and the flux due to the mobility of particles 



subjected to an external force. When there are no exter- 
nal forces except gravity, which also creates the pressure 
gradient, the diffusion flux of a dilute solution reads @ 
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where X is the molar fraction (concentration) of the so- 
lute in the solution, D is the molecular diffusion coef- 
ficient of the solute molecules in the solvent, ar is the 
thermodiffusion constant (ar/T is the Soret coefficient), 
and R = 8.31JK _1 mol~ 1 is the universal gas constant. 
M = M g - AiAf host , where M s and M host are the mo- 
lar masses of solute and solvent molecules, and N\ is the 
number of solvent molecules in the volume occupied by 
one solute molecule in the solution. For specific gases we 
derive Ni from experimental data in Appendix [X] 

When the liquid is saturated with gas bubbles and the 
bubbles are in local thermodynamic equilibrium with the 
solution, the concentration of the solute in the solvent 
equals solubility, X = X^- a \ everywhere in the liquid. 
The thermodynamic equilibrium implies the equality of 
the chemical potentials of the gas dissolved in the liq- 
uid (Mho) an d that of the vapor phase (/if ap ), that is, 
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According to the thermodynamic laws, the 



chemical potential of the vapor phase depends exclusively 
on temperature T and pressure P, //| ap = ^ ap (T, P). 
Hence the solute concentration X is not a free variable, 
but a function of the local temperature and pressure, 
X(f) = (T(r),P(r)); the same holds true for the 
solute flux J(r). In Appendix [Bl of this paper we pro- 
vide the calculation of high-pressure aqueous solubility 
of gases based on a scaled particle theory [l9[ amended 
with implementation of van der Waals' equation of state 
for the vapor phase. 

In the present study we consider the problem of diffu- 
sion transport in a bubbly liquid in the context of real 
geological systems, where pressure varies from one to a 
few hundred atmospheres (this refers to a few kilometers 
deep water column, see Fig. [TJ. To illustrate the ap- 
proach, we will discuss first a more simple case of mod- 
erate pressures P < 100 atm, which allows an explicit 
analytical treatment. The results for the general case, 
obtained numerically, will be also shown. 

According to Eq. (|B12|) in Appendix [B] for moderate 
pressures and far from the solvent boiling temperature 
the solubility depends on T and P as follows 



A( >(T,P)^A(°)(T ,P )^exp 



1 



where To and Pq are reference values, the choice of 
which is guided merely by the convenience reason, and 
X^(Tq,Pq) is the solubility at the reference tempera- 
ture and pressure; the parameter q = —Gi/k-Q is given in 
Tab. |U Pressure is moderate in a sense that one can ne- 
glect van der Waals' effects and the influence of pressure 
on the Gibbs free energy of the cavity formation for the 
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guest molecule in the solvent. The latter is represented 
by the term (-BP/T) in the exponential of Eq. (|BT2]) . 
where constant B ss 10~ 6 K/Pa (see Tab. HJ, and be- 
comes significant only for pressures as high as several 
hundreds atmospheres. The limiting formula is free of 
the molar fraction Y [Eq. (|B13p ] of the solvent molecules 
in gas bubbles because it is assumed to be zero. This as- 
sumption is valid several tens of Kelvins below the water 
boiling temperature, while the systems where liquids are 
at near-boiling conditions, say geysers, are beyond the 
scope of the present study. 

Using VX< ' = l(°)[p- 1 VP-(l + g/T)r 1 VT], one 
obtains for the total flux: 



J 



VP 
~P~ 



— II — CtT 



1) Yl 

TJ T 



Mg 
RT 



(2) 



Note that the flux does not depend on the solute concen- 
tration or its gradient as independent fields. 

Generally, the flux ([2]) possesses a nonzero divergency 
V • </, which implies the existence of sources and sinks 
for the solute mass. Obviously, the gas bubbles serve as 
the mass reservoir, which provides the respective sources 
and sinks. Hence the mass balance reads 



dX b 
dt 



V • J = 0. 



(3) 



where X b is the molar fraction of bubbles in the bubbly 
fluid, that is, in the system comprised by bubbles and liq- 
uid. For the space scales addressed here we do not need to 
consider the microsco pic processes of the bubbles growth 
and clustering as in [2(1 , and the single characteristic, 
-Xb(r), suffices in our case. 

Hence for moderate pressures we obtain the following 
equation for the evolution of a nondissolved gas phase: 



dX h 
dt 
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Eqs. © and (@| are written for the case when the bub- 
bles occupy a small fraction of the fluid volume, which is 
typical for most geological systems, e.g. |a,|6|. Moreover, 
the corrections needed to account for the finite fraction 
of the bubble volume are always quantitative and never 
qualitative. Indeed, the direction of the solute flux is 
merely determined by the factor in the square brackets 
in the right-hand side of Eq. (UJ) , which does not depend 
on the bubble volume. 



III. THE DIFFUSIVE ACCUMULATION AND 
NON-FICKIAN FLUX OF SOLUTES 

In what follows we consider the one-dimensional dif- 
fusion of solutes in bubbly liquids in the presence of 
temperature and density gradients, using for the latter 
quantities some typical geological values. The model of 
one-dimensional diffusion is motivated by the fact that 
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FIG. 2: (Color online) Diagram of diffusive regimes on the 
plane j3—H (/3 quantifies the non-Fickian drift strength) for 
moderate pressures. The critical heights Hi and H2 are given 
by Eqs. ((6)1 and (0, respectively. In the plug-in plots the 
dependence of the solute flux J (blue solid lines) and of the 
bubble growth rate Xt, (red dashed lines) on the depth z below 
the seafloor is shown. 



real geological systems are much more uniform along two 
directions (say horizontal) than along the third direction 
(say vertical). Then we have a laterally uniform sys- 
tem with vertical gradients of its thermodynamic prop- 
erties. For the sake of definiteness we analyze the case 
of seabed sediments, Fig. [TJ where the depth below the 
seafloor is measured by the z coordinate. Then we have 
z-dependent hydrostatic pressure and z-dependent tem- 
perature due to the geothermal temperature gradient: 



P(z) = P + p liqff (z + P), T(z) = T si 



Gz 



Here Pq is the atmospheric pressure, H is the height of 
the pure-water layer above the bubble-bearing porous 
medium, P s f is the temperature of the seafloor and G 
is the geothermal gradient. Nonlinearity of the temper- 
ature profile [2l[ is neglected in our study because it is 
system specific and not a principal ingredient for the phe- 
nomenon we consider. The assumption of a linear tem- 
perature profile is typical for studies on physical processes 
in marine sediments [H, [f| . 

Using the above expressions for the hydrostatic pres- 
sure and geothermal gradient, we obtain for the diffusion 
flux ([2]) for moderate pressures: 



J ~ M< ' 
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FIG. 3: (Color online) Direction of the mass flux and the position of the flux inversion point as the function of the height of 
pure-water column above the sediments. Different curves correspond to different values of the non-Fickian parameter ft for 
methane. Calculations have been performed for the seafloor temperature T s f = 277K and geothermal gradient G = 40K/km (a) 
and G = 60K/km (b) for the general case, without the approximation of moderate pressures. The flux inversion points 
with negative divergency (accumulating the solute gas) are plotted with the solid curves, the points with positive divergency 
( "repelling" the solute gas) are shown with the dashed curves. 
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-CUT 



Mg 
RG' 



Eq. ([5]) with (3 — corresponds to the Fickian diffu- 
sion, J = —DWX, while the coefficient /3 quantifies the 
strength of the non-Fickian contribution to diffusion flux. 
Let us consider the consequences of the above model 
for the typical geological systems, like marine sediments. 
Fig. [2] illustrates the variety of possible diffusive regimes. 
For the gases with positive (3, the gas leaves the sedi- 
ments. In shallow seas the gas diffuses upwards, to the 
sea, in the thin upper layer of sediments and downwards 
below this layer [regime (A) in Fig. [2]. When the sea 
depth H exceeds a certain value Hi, the upper layer 
disappears and the gas diffuses downwards all over the 



seabed [regime (B)]; for moderate pressures the thresh- 
old depth reads, 



H, = 
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Pa 



G (1 + p + q/T s{ ) p liq 
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For gases with negative /3, the direction of the flux very 
deep below the seafloor turns upwards. In deep seas, 
H > Hi [regime (D)], the solute accumulation zone ap- 
pears: The solute gas from the upper and deep layers of 
sediments diffuses to this zone. For shallow seas, H < Hi 
[regime (C)], the solute gas accumulation zone loses its 
contact with the seafloor. In this case the solute gas from 
the region just beneath the seafloor migrates upwards, to 
the sea. 

Thus we reveal a novel and surprising phenomenon — 
the formation of a gas accumulation zone: Instead of the 
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FIG. 4: (Color online) Direction of the mass flux and the 
position of the flux inversion point for carbon dioxide. Cal- 
culations have been performed for the seafloor temperature 
T s f = 277 K and geothermal gradient G = 40K/km for the 
general case, that is, without the the approximation of mod- 
erate pressures. For details see the caption for Fig.[3j 



tendency of the gas to spread uniformly over the system, 
as expected for the common diffusion, the gas concen- 
trates in narrow zones, due to the non-Fickian diffusion 
against the direction dictated by the concentration gra- 
dient. 

Finally, if the positive thermodiffusion is strong 
enough, one observes the regime (E), where no gas- 
retracting zone exists and all the solute migrates up- 
wards, to the sea. Nevertheless, the inhomogeneity of 
the diffusion flux results in the bubble growth zone. For 
moderate pressures, the depth H of the water body be- 
neath which this diffusion regime occurs is bounded from 
above by 



Ho 
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For high pressure, the solubility and diffusion flux are 
more gas-specific and each particular case requires a sep- 
arate analysis. In Figs. [3J and |H the results for methane 
and carbon dioxide are shown for the general case, where 
Eqs. (HJ, (|B12|) . and (|BI3j) have been used without the 
simplifying approximations adopted for the moderate 
pressures. The range of parameters used in Figs. [3J and 
[4] is chosen to be relevant for practical applications: The 
range of pressures for methane corresponds to the possi- 
ble location of the gas-hydrate deposits (methane is one 
of the components of the hydrate) , while the pressures for 



carbon dioxide may be essential for the process of burial 
of industrial CO2 in deep aquifers. 

For methane (Fig. [3]), we consider the upper I — km 
layer of sediments for sea depths ranging from km 
to 5 km with the seafloor temperature T s f = 277 K 
(which corresponds to maximal water density at at- 
mospheric pressure). For a small geothermal gradient, 
G = 20K/km, methane diffuses upwards, to the sea, for 
any value of the non-Fickian parameter 8 and any sea 
depth H . For a larger G, the transport significantly de- 
pends on 8 and H (Fig. [3]). For the typical value of 
G = 40K/km d, § methane diffuses upwards in shallow 
seas and downwards for deep ones; the sea depth H for 
which the flux at z — reverses increases as 8 decreases. 
Remarkably, for the range of the non-Fickian parameter 
0-3 < 8 < 0.9, the gas accumulation zone (which blocks 
the release of methane into the sea) is located in the 
seabed at the depth ranging from 0.7 km to 1.6 km (solid 
lines in Fig. [3^). For 8 < — 1, this zone appears in sedi- 
ments under deep water bodies, H > 3.9 km. In Fig. [3j3, 
where G = 60K/km, one can see how the diffusion 
regimes are affected by further increase of the geother- 
mal gradient. We wish to stress, that the reported phe- 
nomenon originates exclusively due to the non-Fickian 
diffusion. The oversimplified models based on the Fick- 
ian law with 8 = 0, e.g. [!, Q, cannot describe the novel 
effect of the formation of gas accumulation zone. 

The main difficulty of the practical application of our 
theory is the lack of precise data for the thermodiffu- 
sion constant, that is, the lack of a trustworthy value 
for the coefficient 8. For instance, the authors are not 
aware of any experimental data on the thermodiffusion 
of methane in water. Theoretical studies, e.g. |22j, do 
provide the first-principle expressions for computing the 
thermodiffusion constant. These, however, require the 
knowledge of the intermolecular potentials and structural 
properties of the solvent. Hence, to obtain a reliable 
estimate for the non-Fickian parameter, extensive nu- 
merical simulations by means of Molecular Dynamics or 
Monte Carlo are needed, which is beyond the scope of the 
present study. Presently, however, we can prospect the 
value of «t from a simplified model of a large Brownian 
chemically-inert particle in non-isothermal liquids. With 
this model, cut is controlled merely by the volume per 
one guest molecule in the solution and the guest molecule 
mass; employing Buckingham's 7r-theorem of dimensional 
analysis [23J one can make an interpolation of the data for 
dilute aqueous solutions of paraffin alcohols: olj- ~ 1.5 for 
methanol 24 1, ax ~ 3.0 for ethanol [25] . and ar ~ 4.5 
for isopropanol 26]. This interpolation yields ay ~ 1.8 
for methane. Thus, /3 CH4 w —2.5 for aqueous solutions in 
the presence of geothermal gradient G = 40K/km, and 
^CH* ~ _ 2 .3 f or G = 60K/km. 

For the case of carbon dioxide the non-Fickian drift af- 
fects the solute flux much more weakly than for methane 
(Fig. 2]). Under shallow water bodies, CO2 diffuses up- 
wards in upper layers of sediments, and downwards for 
deep layers. 
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IV. CONCLUSION 

We have developed the theoretical description of the 
process of diffusive migration of a dissolved gas in bub- 
bly liquids, where bubbles are immovably trapped by a 
porous matrix, as in a seabed or terrestrial aquifers. The 
effect of temperature inhomogeneity across the system 
and the gravity force are accounted for. The theory is 
employed for treatment of typical geological systems with 
the hydrostatic pressure distribution and the geothermal 
gradient. 

For non-high pressure the diffusion flux and bubble 
mass evolution are governed by Eqs. © and ((3]), respec- 
tively; Fig. [5] presents the diagram of diffusive regimes 
derived for these equations. For high pressure the sys- 
tem is comprehensively described by Eqs. (JXJ) , (|B12|) . and 
(IE>13|) with parameters provided in Tab. U for four typical 
gases: nitrogen, oxygen, carbon dioxide, and methane. 
The high-pressure solubility (|B1 21) has been derived from 
the standard scaled particle theory pl| with adoption of 
van der Waals equation of state. 

The "non-Fickian" corrections — thermodiffusion and 
gravitational segregation — appear to play an essential 
role in migration of methane in the seabed, even being 
able to create the gaseous methane accumulation zone in 
sediments (Fig. [3]). For carbon dioxide the non-Fickian 
effects are not as significant: they do not cause the quali- 
tative change of behavior but can only shift the boundary 
of the CC>2-capture zone by not more than 75 m, which 
is 15% of the characteristic depth (Fig. [4j. Unfortu- 
nately, the precise values of the thermodiffusion constant 
of aqueous solutions of methane or carbon dioxide are not 
found in the literature, and one can only estimate their 
values as we do in this paper. The practical significance 
of the thermodiffusion effect, which has become appar- 
ent for the evolution of geological systems, necessitates 
the experimental determination of the thermodiffusion 
constant for aqueous solutions of methane and carbon 
dioxide. 
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Appendix A: Calculation of Ni and M from 
experimental data 

The ratio of the volume occupied by one mole of solute 
molecules in the solvent, vf iq , and the molar volume of 
this solvent, Uii q , [which is Ni = ( v f iq /vu q ) in SecQl] can 
be precisely derived for CO2-H2O and CH4-H2O sys- 
tems from the dependence of the solution density on its 



concentration, which is available in the literature |27H29| . 
The molar volume of solution is u so i = v n q (l — X) + vf iq X, 
and the density is 
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host 
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(1 - X) + M S X 



"sol 



Psol(* = 0) 



X + 0(X 2 ) 



where M host and M s are the molar masses of the solvent 
and solute molecules, respectively; p so \(X — 0) = pn q is 
the density of a pure liquid (solvent). Hence, 



vn q M host 
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sol 



PHq 
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Pliq 



sol 
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The results of the implementation of the latter relation 
to the experimental data from [271 . l28j | are summarized 
in the three last columns of Table Q] 



Appendix B: Aqueous solubility of certain gases 

Here we wish to present some detail of our calculations, 
outlined in previous sections. To evaluate the solubility 
at high pressure we make a straightforward generalization 
of the standard scaled particle theory. Namely, instead of 
the ideal-gas equation of state we use the van der Waals 
equation for the vapor phase of the dissolved gas. We find 
the chemical potential of the solute and calculate solu- 
bility for the gases, most important for applications — 
oxygen, nitrogen, carbon dioxide, and methane. 



1. Chemical potential in the vapor phase 

To introduce notations and recall basic relations we 
start from an ideal gas. The Helmholtz free energy, 
F(T,V,N), the Gibbs free energy, G(T,P,N) = F + 
PV = F - V(dF/dV), and the chemical potential, 
//f ap (T,P) = dG(T,P,N)/dN, read for an ideal gas, 



F = 

G = 

u g = 

r^vap 
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(Bl) 
(B2) 
(B3) 



where N is the number of molecules, V is the vol- 
ume, P is pressure, k B is the Boltzmann constant, A = 
yj2Trh 2 /mk B T is the thermal de Broglie wavelength of a 
gas molecule of mass to, h is the Planck constant divided 
by 2n and Z(T) is the partition function of the internal 
degrees of freedom of the molecule. It depends on the 
molecular structure and temperature. 
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Now, we address the chemical potential of real gases, 
taking into account the finite volume fraction of the gas 
molecules and their attractive interactions. The simplest 
way to do this is to use the van der Waals model. We as- 
sume that the heat capacity cy does not depend on tem- 
perature. This is justified for the range of temperatures 
of interest for many gases, including methane. Then the 
Helmholtz and Gibbs free energies read, e.g. [30(, 

F(T,V,N) = F id .-Nk B Tln(l-nb)-Nna,(B4) 
G(T, P, N) = G id . - Nh B T In (1 - nb) 



- 2Nna + Nk B T 



nb 



1 — nb 



(B5) 



Here -Fid. and Gid. are the respective ideal parts and a 
and b are the van der Waals constants (see Tab. U for the 
specific values). Note that above the critical point of a 
gas (T c = — 82.7°C for methane), n = n(P) is a univocal 
function of pressure, which may be obtained from the van 
der Waals equation of state: 



P 



nk B T 
1-nb 



an 



(B6) 



Finally, the chemical potential of the van der Waals gas 
reads, 



t ^ KP (T,P) = k B T 



In 



nA 3 

W) 



- ln(l - bn) 
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This result can be found in 



2. Solubility of gas: Scaled particle theory 

To evaluate the solubility of a gas in a solvent (water) 
we employ scaled particle theory (see e.g. the review [19[ 
for detail). When the concentration of the dissolved gas 
reaches solubility, the system is in thermodynamic equi- 
librium, so that the chemical potential of the gas in the 
vapor phase ^| ap is equal to this in the solution /xf iq . Ac- 
cording to the scaled particle theory [19| , 

/xf iq = G C + G l + k B Tln{A 3 /Z) + k B Tln(X N A /v liq ), 

(B8) 

where A and Z are the de Broglie wavelength and the par- 
tition functions per molecule for the internal degrees of 
freedom for the solute, X is the molar fraction of the gas, 
G c is the work of creation of a cavity for a guest molecule 
in the solvent, Gi is the interaction energy between a so- 
lute molecule and the surrounding solvent molecules, Aa 
is the Avogadro number, and vy lq is the molar volume 
of the solvent particles (e.g. «h 2 o = 17.95 cm 3 /mol at 
atmospheric pressure and T = 300 K) . The cavity for- 



mation work reads [1 

G r 



k B T 



-ln(l - 


-y) + j 






}-y 




yP 





'Ay cr g 



Cliq 



n\i q k B T Vciiq, 



A(y)+B(y)-, (B9) 



where cri iq and a g are respectively the effective diameter 
of the solvent (liquid) and solute (gas) molecules, ni; q is 
the number density of solvent, and y = (7r/6)nii q of iq (e.g., 
y B2 o = 0.371 at T = 300 K 19]). There exist numerous 
models for the interaction energy Gi (e.g. [l9|, HH); for 
the present study it is enough, however, to mention that 
it is almost independent of pressure and temperature. 

For the vapor phase, the chemical potential is given 
either by Eq. (|B3[) for an ideal gas or by Eq. (|B7|) with 
(|B6|) for n = n(P). 



/'v 



k B T\n(A 3 /Z) + k B Tln(f/k B T), 



(B10) 



where / is the fugacity of the gas molecules in the vapor 
phase. It is equal to the pressure for a one-component 
gas, and to the partial pressure for a gas mixture. In 
the case of interest, / = P(l — Y), where Y is the molar 
fraction of solvent vapor in the gaseous phase; it is de- 
termined by Eq. (|B13[) below. Equating Eqs. (|B8|) and 
(|B10|) for the chemical potential, one finds the equilib- 
rium molar fraction of the dissolved gas, that is the solu- 
bility, For an ideal gas one obtains, using Eq. (IB10|) 
together with the ideal gas equation of state 



X (o) 



. (l-Y)Pvy lq 

RT 
(1 -Y)P v liq 
RT 



'■ exp 
exp 



G c + Gi 
k B T 



-A{y) — B— - 



(Bll) 



For a van der Waals gas, the real gas equation of state, 
Eq. (|B6[) . is to be employed. This yields 



X (o) 



■ exp 



(1 - Y)nv Uq 
(1 - nb)N A 

_ (1 - Y)nv Vlq 
(l-nb)N A 



G c + G t 2an 



exp 



k B T 

My)- 

Ian 



k B T 
P 



B- 



1 — nb 
G, 



T k B T 
1 



k B T 1 



no 



(B12) 



When liquid with a dissolved gas is in equilibrium with 
its vapor phase, the vapor, including the vapor in bub- 
bles, contains both components. Let us evaluate the 
solvent molar fraction Y in the vapor phase. The en- 
thalpy variation dH = c P dT + VdP = T dS + VdP. 
Hence, H = H + cp^) dT x + V(T, P) dP, S = 

So + f^cpiT^T^dTi, and Gibbs free energy G = 
H~TS = H + J^cp(T 1 )dT 1 +Jp g V(T,P)dP-TS - 
Tflcp^T^dT,. 
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TABLE I: Parameters of gases and their aqueous solutions. Values of a, B, and Gi are derived from the theoretical fitting of 
the experimental data presented in Fig. [5] Parameters in the three last columns are derived from experimental data [27l. [28l] 
with Eq. (TAT]) . 





m 6 Pa 
mol 2 


k m 3 
b, W- 3 — 
mol 


0, A 


My) 


B, 10-^ 


ks ' 


1 dp so i 
po dX 


vf- 

Ni = -is 


M, -i- 
mol 


2 


0.1378 


3.183 


3.03 


6.38 


1.056 


-873 








N 2 


0.1408 


3.913 


3.15 


6.77 


1.182 


-781 








C0 2 


0.3640 


4.267 


3.02 


6.36 


1.050 


-1850 


0.54 


1.91 


9.7 


CH4 


0.2283 


4.278 


3.27 


7.18 


1.321 


-1138 


-1.35 


2.24 


-24.3 


H 2 






2.77 















For liquid water, cp(T) is nearly constant and 
G!i° q st = fl£$ + c h P f q (T - To) + ^(P - Po) 



iiq,o 



T5, 



host 
liq.O 



for the vapor phase, 



host 
-P.liq 



T 

Tin — 



T n 



where 



„H a O 
-*P,liq 



- r H 2 
c P,vap 

9.09, the 



„h 2 o 



.H 2 



'P,Uq> -P,vap/ fc B = 4.09, 

enthalpy of vaporization 

AH^ 2 °/k B = 4892K at T = 373K and P = latm. 
Eq. (|B13p provides the approximate value of the sat- 
urated vapor pressure which is nearly indistinguishable 
from the ex per imental value and the known empiric for- 
mulae (e.g. |3a|). 



y^fhost rrhost 

^vap vap,0 



c*% p (T -T ) +k B T In ^~ 



YP 
Po 



rp Qhost 

± °vap,0 



T 

c p"vap^ m 77-T 



„host 



Here Y appears owing to the fact that the part of en- 
tropy related to the translational degrees of freedom is 
k B ln(V/FA 3 ) but not k B ln(V^/A 3 ) as for a pure steam. 

In equilibrium, the chemical potential in the liquid and 
vapor phases are equal, G^f = G[|° st , and 



Y 



h 2 o 





fcs 




exp 



vn q (P - Po) 



AH, 



h 2 o 



- Ac P 2 °T 



RT 
1 1 



(B13) 



Table Q] presents parameters a, B, and Gi derived 
from experimental data (31rl33 | with the employment of 
Eqs. (|B11I) and (IB13I) for an ideal gas and Eqs. (IB12I) and 
(|B13|) for van der Waals' model. The agreement between 
the theory and experimental data can be assessed from 
Fig. [5j Notably, van der Waals' model provides an accu- 
rate description in the entire range of parameters relevant 
to our consideration. One should keep in mind that van 
der Waals' model provides only a qualitative description 
for liquid-gas transition and, therefore, the formulae we 
present are not applicable for the description of dissolu- 
tion of any gas at its liquefaction conditions. 
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FIG. 5: (Color online) Solubility of O2 (a), N2 (b), CO2 (c), and CH4 (d) in water at specified pressure. Squares represent 
experimental data f [3ll. I32T] for nitrogen, oxygen, and carbon dioxide and [33, [34[ for methane); red dashed lines: results of the 
scaled particle theory for vapor phase assumed to be an ideal gas, Eq. (|Blf [) ; blue solid lines: the scaled particle theory with 
van der Waals' equation of state for the vapor phase, Eq. (|B12|) . 
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